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ABSTRACT 

We present a new solution of the direct problem of planet transits based on transfor- 
mation of double integrals to single ones. On the basis of our direct problem solution 
we created the code TAC-MAKER for rapid and interactive calculation of synthetic 
planet transits by numerical computations of the integrals. The validation of our ap- 
proach was made by comparison with the results of the wide-spread Mandel & Agol 
(2002) method for the cases of linear, quadratic and squared root limb-darkening laws 
and various combinations of model parameters. For the first time our approach allows 
the use of arbitrary limb-darkening law of the host star. This advantage together with 
the practically arbitrary precision of the calculations make the code a valuable tool 
that faces the challenges of the continuously increasing photometric precision of the 
ground-based and space observations. 

Key words: methods:analytical - methods:numerical - planetary system - bina- 
ries: eclipsing 



1 INTRODUCTION 

The determination of geometric parameters of extrasolar 
planets has an essential role in inferring their densities and 
hence their compositions, masses and ages. This information 
leads to refinements of the models of planetary systems and 



yields important constraints on planet formation (Cody & 



Sassel ov 2002; Hubbard et al.|2001[|Seager fc Mallen-Ornelas 



2003) 



During the last several years there is a sharp rise in the 
detections of transiting extra-solar planets (TEP) mainly by 
the wide-field photometric variability surveys: (i) ground- 



based observations as Super WASP (Pollacco et al 



HATNet (Bakos et al. 20041, OGLE-III (Udalski et al 



2006), 



2002). 



TrES (|Alonso et al.|2004| 7etc; (ii) space missions as CoRo T 



(Bagl in et al.|2006| l and Kepler ( |Borucki et al.|2010a l. 

The Kepler mission produced a real bump of the num- 
ber of exoplanet candidates, dozens of them yet confirmed 



| Koch et al.|2010||Borucki et al.|20l"0b"l|Dunham et al.|20f0 



Latham et al. 2010 etc.). The space-based missions, espe- 



cially Kepler, have the photometric ability to detect even 
transiting terrestrial-size planets ( Steffen et al.|20T2 l. 

The recent increasing number of the planet-candidate 
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discoveries and the increasing precision of the observations 
allow to investigate fine effects such as: 

(a) rotational and orbital synchronization and alignment 
( |Winn et al.|[20TT] [Hibrard et al.||2lH0] ); 

(b) zonal flows and violent atmospheric dynamics due to 
large temperature contrast between day-sides and night- 
sides of the planets ( Burrows et al.|2010 |; 

(c) departures from sphericity of the planets and spin pre- 
cession ( |Carter fc Winn||2010[); 



(d) rings and satellites ( |Hui fc Seager||2002 
Schneider|2006l ); 



Arnold & 



(e) stellar spots (Rabus et al.|2009 Dittmann et al.|2009 



|Hebrard et al.|2010| ) 

(f) dependence of the derived planetary radius from the 



limb-darkening coefficients ( Kipping fc Bakos|2011[); 

(g) irradiation by the parent star (Seager & Whitney 



20001, planet transmission spectra (Seager & Sasselov|2000 1, 



atmospheric lensing due to atmospheric refraction (Hui & 
Seager 2002), Rayleigh scattering, cloud scattering, refrac- 
tion and molecular absorption of starlight in the planet at- 
mosphere ( Hubbard et al.|2001 1, etc. 



The study of such fine effects requires very precise meth- 
ods for determination of parameters of the planetary systems 
from the observational data. 

Recently, we established that the synthetic light curves 
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generated by the widely-used codes for planet transits de- 
viated from the expected smooth shape. This motivated us 
to search for a new approach for the direct problem solu- 
tion of the planet transits. We managed to realize this idea 
successfully for the case of orbital inclination i = 90° and 
linear limb-darkening law ( Kjurkchieva fc Dimitrov||2012[ ). 
This paper presents continuation of the new approach for 
arbitrary orbital inclinations and arbitrary limb-darkening 
laws. 



2 METHODS FOR SOLUTION OF THE 
PLANET TRANSIT PROBLEM 

2.1 Previous approaches 

The global parameters of the configurations of TEPs are 
different from those of eclipsing binary systems. The main 
geometric difference is that the radii of the two components 
of TEPs are very different, while those of EBs are compa- 
rable. As a result, almost all models based on numerical 
integration over the stellar surfaces of the components give 
numerical errors, especially around the transit center. The 
overcoming of this problem required specific approaches and 
models for the study of TEPs. 

The first solution of the direct problem of the planet 



transits is that of Mandel & Agol (2002). They derived an- 



alytical formulae describing the light decreasing due to cov- 
ering of stellar disk by a dark (opaque) planet in cases of 
quadratic and nonlinear limb-darkening laws. The formulae 
oflMandel & Agoll (120021 further M&A solution) contain 



several types of special functions (beta function, Appel's hy- 
pergeometric function and Gauss hypergeometric function, 
complete elliptic integral of the third kind). To generate syn- 
thetic transits [Mandel fc AgoI| ( |2002[ ) created IDL and For- 
tran codes occultsmall (for small planet), occultquad 
(for quadratic limb-darkening law) and OCCULTNL (for non- 
linear limb- darkening laws) that are based on numerical cal- 
culations of the special function values. These codes are 
widely used by many investigators for analysis of observed 
transits and improved later by different authors. 

In the meantime ISeager & Mallen-Ornelas ( 2003 1 ob- 



tained analytical solution for the particular case of total 
transit and uniform stellar disk (i.e. neglecting the limb- 



darkening effect), Knutson et al. (20071 made calculation 



of the secondary planet eclipses while Kipping (2008 20101 



studied the problem of eccentric orbits. 

The second direct problem solution for the planet tran- 



sits was made by Gimenez (2006). He derived analytical for- 



mulae for the computation the light curves of planet' transits 
for arbitrary limb-darkening laws. This approach is similar 
to that of the Kopal's a n functions and the derived formu- 
lae contain different special functions (elliptic integrals of 
the first, second and third order). 

Most of the codes for inverse problem solutions of planet 
transits are based on the |MfcA| soluti on. For instance, the re- 



cent packages TAP and autoKei^] ( G azak et al.|2012 1 



Pal 



et al. (20101 estimated the fit quality of the inverse problem 



corresponding to the |Mfc A~] solution. It should be noted that 
the inverse problem solution of the planet transits is not a 



trivial task. The known codes for stellar eclipses are not ap- 
plicable for the analysis of most planet transits due to the 
non-effective convergence of the differential corrections in 
cases of observational precisions poorer than 1/10 the depth 
of pla net transit. EBOP ( [Etzel|1975[ |1981| |Popper fc Etzel 



1981 1 is the only model for EBs, which heavily modified ver- 



sion JKTEBOP ( [Southworth|2008||Southworth, Maxted, fc 
Smalley|2004a|b I can be applied successfully to TEPs 



Recently, solutions of the whole inverse problem based 
on simultaneous modeling of photometric and spectral data 
of exoplanets performed using the Markov-chain Monte 



Carlo (MCMC) code were proposed (Collier Cameron et 
al. 2007 Pollacco et al. 2008 etc.). But their subroutines 



for fitting of the transits also use the |MfcA| solution for 
quadratic limb-darkening law. For instance, the subroutine 



exofast_occultquad of the package exofast ( Eastman et 
al.|2013 1 is a new improved version of OCCULTQUAD 



An opportunity to generate synthetic transit light curve 
as well as to search for fit to own observational data is pro- 
vided from the website Exoplanet Transit Database]^] (ETD). 
It is based on the OCCULTSMALL routine of the IMfcAl solu- 
tion that uses the simplification of the planet trajectory as a 



straight line over the stellar disk (Poddany, Brat, & Pejcha 



2010) 



2.2 The new solution of the direct problem: how 
to use the symmetry of the problem to reduce 
the double integral to a single one 

Let's consider configuration from a spherical planet with 
radius R p orbiting a spherical star with radius R s on circle 
orbit with radius a, period P and initial epoch To. Let's the 
line-of-sight is inclined at an angle i to the orbital plane of 
the planet. 

Usually, it is assumed that the limb-darkening of the 
main-sequence stars may be represented by the linear func- 
tion 



/(/x) = / [l -«(!-/*)] 



(1) 



where Jo is the light intensity at the center of the stellar disk 
depending on the stellar temperature, fj, — cos 9 and 8 is the 
angle between the normal to the current point of the stellar 
surface and the line of sight. 



Claret ( 2000 I found that the more accurate limb- 
darkening functions are the quadratic law ( |Kopal|1950[ ) 

/(/x) = J„ [1 - ui(l - (J.) - ua(l - m) 2 )] (2) 



and 



J(aO = Jo 



nonlinear" law 

4 

1 



3=1 



(3) 



which is a Taylor series in fj, to fourth order in 1/2. 

Square-root law is proposed by Diaz-Cordovcs & 



Gimenez ( 1992 1 



J(p) = I [l-ui(l-A»)-Ua(l-VM)] (4) 

while Klinglcsmith & Sobieski ( |1970[ ) proposed logarithmic 
law for early stars 



http:/ /ifa. hawaii.cdu/uscrs/zgazak/IfA/TAP.html 



http: / /var2. astro.cz/ETD/protocol.php 
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= J [1 — ui(l — — ua/iln/i] . 



(5) 



Our solution can be applied for arbitrary limb- 
darkening law 



(6) 



where /(itj,/i) is an arbitrary function of /x. 

The possibility to use an arbitrary limb-darkening law 
is one of the main advantages of our approach. 

The luminosity of the planetary system at phase <p out- 
of-transit is 



L(ip) = L s + Lp(ip) 



(7) 



where L s is the stellar luminosity and L p (<p) is the planet 
luminosity. The phase tp is calculated by the period P and 
the initial epoch HJD(mm). 

The planet luminosity L p (p) is variable out of the 
eclipse (because more of the day-side of the planet is vis- 
ible after the occultation while more of the night-side of the 
planet is visible before the transit). Due to the relatively 
small size and low temperature of the planet, it can be as- 
sumed that its disk is uniform and its luminosity does not 
change during the transit, i.e. 



L p = irRpIp 



where I p depends on the planet temperature. 
The luminosity during the transit is 

L(ip) = L s + L p — J(ip) 



(8) 



(9) 



where J(<p) is the light decrease due to the covering of the 
star by the planet. It might be expressed in the form 



J(<P) = 



(10) 



where the integration is on the stellar area S oc (<p) covered by 
the planet. Hence, the solution of the direct problem for the 
planet transit is reduced to a calculation of a surface integral. 
It is not a trivial task due to the nonuniform-illuminated 
stellar disk. 

If we assume the out-of-transit flux to be F ou t — l, then 
its decreasing during the transit is described by the expres- 
sion 



L s + Lp - J(ip) 



(11) 



L s + Lp 

For the next considerations we use coordinate system 
whose origin coincides with the stellar center. The axis z is 
along the line-of-sight and the xy plane coincides with the 
visible plane. We choose the axis y to be along the projection 
of the normal to the orbit on the visible plane. 

Taking into account that the stellar isolines of equal 
light intensity are concentric circles with radius r we may 
calculate the light decrease as a sum (integral) of the contri- 
butions of differential uniformly-illuminated arcs with cen- 
tral angles 2j T and area ds — 2j T rdr (Fig. W. In this way 
we transform the surface integral ( 10 1 to a linear one 



1-maxO) 
r m i„(¥>) 

where 



(12) 
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Figure 1. The visible plane: geometry of the partial transit (Case 
A.l) 



/'• : 



cos - 



(- 



(13) 



The integrand of our main equation ( 12 1 is different 
from that of the main equation (2) of the M&A solution. 
As a result, the methods of numerical calculation of these 
integrals are different. 

The integration limits r m i n (v?) and r max ((p) are the ex- 
tremal radii of the stellar isolines that are covered by the 
planet at orbital phase ip. These limits depend on the con- 
figuration parameters. 

It is appropriate to assume the separation a as a size 
unit and to work with dimensionless quantities: relative ra- 
dius of the planet r p = R p /a; relative radius of the star 
r s — Rs/a; relative radius of the stellar isoline q — r/a. 

Then the equation for the light decrease during the tran- 
sit can be rewritten in the following form: 



F{<p) = F ou 



J{fp) 



■Kkr% + L s /[I a 2 ] 
where 

9maxO) 

J(tp)= J f{u- s ,^)2^ q {ip)qdq 

tmin(f) 

k = / p //o = [Tp/Tof 



(14) 



(15) 



(16) 



Assuming the center of the transit to be at phase 0.0 
we consider only the phase interval (0, 0.5) because in the 
case of spherical star and planet the transit light curve in 
the range (-0.5, 0) is symmetric to that in the range (0, 0.5). 

For a circular orbit the coordinates (in units a) of the 
planet center at phase ip are 



xo(ip) — sin(27T<^) 
yo(ip) = cos i cos(2-Kip). 



(17) 



The coordinates xi^i'p), yi^i'p) (in units a) of the in- 
tersection points of the stellar brightness isoline with radius 
q and the planet limb are respectively (Fig. [T]): 



xx(cp) = 
x 2 (p) = 



c 2 (p)x () (p) + y (p) v / 4b 2 (p)q 2 - c 4 (y) 
2b\p) 

c 2 {p)x ((p) - yo((p)y/4b 2 (p)q 2 - c 4 (y) 
2b 2 {p) 
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Figure 2. Geometry of partial transit in Case A. 2 



yi(y) = 

2/2 (<£) = 



c 2 (y%oQ) ~ x (<p)y/4b 2 (<p)q 2 - c 4 (o3) 
2& 2 (v) 

c 2 (<p)y (<p) +x Q (<p)^/4b 2 {<p)q 2 - c 4 (o3) 



2& 2 (o9) 

where we have introduced the designations 



(18) 



b 2 ((p) = xq(>p) + 2/0(09) = 1 — cos 2 (27r<p) sin 2 i 

c 2 (o9) = ? 2 + 6 2 M-r 2 . (19) 

Further we will derive the expressions for the angle 



27 q (<p) and the integration limits in equation ( 15 1 for differ- 
ent combinations of geometric parameters and at different 
phases. 




Figure 3. Geometry of partial transit in Case A. 3 
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2.2.1 Case A: Partial transits 

The transit is partial if the orbital inclination is into the 
range i\ < i < i 2 where 

i\ = arccos(r s + r p ) 

22 = arccos(r s — r p ). (20) 

The phases of outer contacts star-planet are — ooi (be- 
ginning of the transit) and +<pi (end of the transit) where 

1 . ^(r s + r p ) 2 -cosi 2 

03i = — arcsm — . (21) 

Y 2tt sini v ; 

The partial transit occurs into the phase range [0, 09i]. 
For the sake of brevity we will not write further the depen- 
dence of x ,y ,x 1 ,x 2 ,yi,y2,b,c on tp. 

(A.l) If xo — \Ap — 3/0 > r p (Fig. II) the light decreasing 
is calculated by ( 15 1 where limits ana expression for 2^ q (tp) 
are given in Table 1 (Case A.l). 

(A.2) If x - 



/rf-Vo < r p < x Q + ^Jr 2 -y' 2 then the 
integral J(<p) in ( 15 I can be presented as a sum of two in- 
tegrals Ji (09) + J2 (09) (Fig. [2| whose limits and expressions 
for 27q(<p) are given in Table [I] (Case A.2). 

(A. 3) If xo + \/fp — Do < r p the integral J(ip) (Fig. is 
a sum of three integrals Ji(ip) + Jvif) + ^(v) whose limits 
and expressions for 27 q (o9) are given in Table [l] (Case A. 3). 



2.2.2 Case B: Total transits out the stellar center 

If the orbital inclination is into the range i 2 < i < 13 
where 



Figure 4. Geometry of the total transit, Case B 

then the transit develops from partial to total. In this case 
the planet does not cover the star center at any phase (Fig. 

0- 

The phases of the inner contact star-planet are —ip2 
(during planet' entering) and +092 (during planet' exit) 
where 



092 = — arcsm — — : . (23) 

27T sin i 

Into the phase range [052,091] the transit is partial and 
geometry is similar to the Case A. Into the phase range 
[0, 092] the transit is total. 



(B.l) If yo > r p then the integral J((p) is calculated by ( 15 1 
where limits and expression for 27 q (o9) are given in Table]]] 
(Case B.l). 

(B.2) If yo < r p the integral J(<p) is a sum of three integrals 
Ji(ip) + ^2(09) + Js(oo) whose attributes are given in Table [l] 
(Case B.2). 



2.2.3 Case C: Total transit through the stellar center 

If the orbital inclination is into the range is. < i < 90° the 
transit develops from partial to total and finally the planet 
covers the stellar center. 

The phases at which the planet limb touches the stellar 
center are —093 (before the transit center) and +093 (after 
the transit center) where 



iz — arccos (r p ) 



(22) 



¥>3 



2ty 



(24) 
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Table 1. Integration limits and expressions for the integrands used for our direct problem solution 



Integrals 


Cases 
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9max(</>) 
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r s 
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Figure 5. Geometry of the total transit for yo < r p , Case B.2 



Into the phase range [ip2,v>i] the transit is partial and 
the geometry is similar to the Case A. 



If yo > r p then the integral J(ip) is calculated by (151 
where limits and expression for 27 q (<^) are given in Table]]] 
(Case A.l). 

If yo < r p and xo + y/rp — j/q > r s then the integral J(v>) 
is presented as a sum of two integrals Ji(<p) + Jzif) whose 
attributes are the same as those of the Case A.2 (Table 

If yo < r p and xo + \J r p — y\ < r s then the integral 
J{ip) is presented as a sum of three integrals Ji (</?) + J2 (tp) + 
J2,{<p) whose attributes are the same as those of the Case A.3 
(Table [l]). 

Into the phase range [753 , ip2~\ the transit is total, outside 
the stellar center. The geometry is similar to Case B. 



where limits and expression for 27 q (y) are given in Table [l] 
(Case B.l). 

If yo < r p then the integral J(v>) is presented as a sum 
of three integrals Ji (ip) + J2 (<p) + J3 (v) whose attributes are 
the same as those of the Case B.2 (Table [TJ. 

Into the phase range [0, ipz] the planet covers the stellar 
center and there are three subcases depending on the phase 



(pi = — arctan(cosi) 

27T 

that corresponds to the moment when xq — yo- 



(25) 



(C.l) If (£4 < if3 then into the phase range [^4,1^3] (for 
which xo > yo) the integral J(v>) is presented as a sum of 
six integrals Ji((f) + Mv) + Mv) + Mv) + M^v) + M'fi) 
(Fig. |6| whose attributes are given in Table [l] (Case C.l). 

(C.2) If tfi < (f3 then into the phase range [0, ^4] (for 
which xo < yo) the integral J(<p) is presented as a sum of 
six integrals Ji(<p) + Mv) + J'a(v) + Ji(<p) + Jsif) + Mv) 
(Fig. [7| whose attributes are given in Table [l] (Case C.2). 

(C.3) If > 993 then into the whole phase range [0, ^3] is 
fulfilled xo > yo and the integral J(f) is presented as a sum 
of six integrals Ji (ip) + J2 (<p) + J3 (<p) + J4 (<p) + J5 (<p) + Je (<p) 
(Fig. [6| whose attributes are the same as those of Case C.l. 

Note: The condition 954 > ip$ is satisfied for 13 < i < 14 
where 



1-4 



(26) 



If yo > r p then the integral J(ip) is calculated by (151 



Finally, it should be noted that for the arbitrary limb- 
darkening law (JsJ) the stellar luminosity L s is given by the 
integral 
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Figure 6. Case C.l and C.3, total transit through stellar center 
and xq > j/o 




Figure 7. Case C.2, total transit through stellar center and 

xo < yo 



hf{u- 3 ,^)2^rdr. 



(27) 



It has analytical solution for all known limb-darkening func- 
tions. Particularly, for the wide-used limb-darkening laws 
(equations from |l]) to |5|) the stellar luminosity is calcu- 
lated by the formulae 



is 



= nR s Jo ( 1 



= 7ri? s Jo I 1 



ttR e Jq I 1 



nR s Jo ( 1 



= 7T R B Jo I 1 
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3 THE CODE TAC-MAKER 



linear law (28) 
quadratic law 
squared-root law 
logarithmic law 
^ ) "jwmliiiciv" law 



The integral in equation |15j) cannot be solved analytically. 
That is why we had to carry out a numerical solution of this 
integral. 

For this reason we wrote the code TAC-maker (Transit 
Analytical Curve) whose input parameters are: 

• radius of the orbit a; 

• Period P and initial epoch To; 

• radius of the star R B ; 

• radius of the planet R p ; 

• orbital inclination i; 

• temperature of the star T s ; 

• temperature of the planet T p ; 

• coefficients of the limb-darkening uf, 

• step in phase PH Btcp ; 

• parameter of precision PPtac of the numerical calcu- 
lations of the integrals. 

The code TAC-maker provides the possibility to 
choose the limb-darkening law from a list of the known 
wide-spread functions (1-5), or to write arbitrary function 
/(«j, fi). This is a significant advantage of the proposed ap- 
proach as even now the accuracy of the most used quadratic 
limb-darkening law is worse than the achieved from Kepler 
for large planets with R p > 0.04Ji s ( jEastman et al.|2013[ ). 

Moreover, the code TAC-maker allows to obtain the 
stellar limb-darkening coefficients from the transit solution 
and to compare them with the theoretical values of |Claret] 
(2004) calculated for different temperatures, surface gravi- 
ties, metal abundances and micro-turbulence velocities. 

The code is written in Python 2.7 language with a 
Graphical User Interface. At the very beginning the code 
checks which condition (Case A, Case B, or Case C) is sat- 
isfied for the given combination of configuration parameters 
(stellar and planet radii and orbital inclination). After that 
the code calculates the characteristic phases tpi for the re- 
spective case. Further the code makes numerical calculation 



of the integral ( 15 \ for the current phase and chosen function 



/(uj,/x) using the SCIPY package. Finally, the code repeats 
the procedure for each phase of the corresponding phase 
ranges. The output results flow as data file (phase tp, flux 
F). 

The code allows to search for solutions of observed tran- 
sits by the method of trials and errors varying the input 
parameters. The data file might be in format magnitude or 
flux and HJD or phase. An estimate of the fit quality is the 
calculated value of \ 2 - Moreover, the two plots of the current 
solution showing the observational data with the synthetic 
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transit as well as the corresponding phase distribution of the 
residuals allow fast finding of a good fit. 



3.1 Validation of TAC-maker 

To validate our approach we used comparison with the wide- 
spread M&A solution, particularly we compared the syn- 
thetic light curves generated by the code TAC-maker and 
those produced by the code OCCULTNL for the same config- 
uration parameters and the same limb-darkening law. For 
this purpose we applied the freely available version of the 
last code without any changes, particularly with its default 
precision, while our code worked with the default precision 
of the SCIPY package for the numerical calculations of the 
integrals. 



3.1.1 Comparison for linear limb- darkening law 

We established that for linear limb-darkening law the two 
synthetic transits coincide (Fig. [8j top). However, the de- 
tailed review reveals the meandering course of the Mfc A| so- 
lution around the smooth course of the TAC-maker transit 
curve (Fig. [81 second panel). The phase derivatives of the 
fluxes (Fig.j8j third panel) exhibit more clearly which one 
of the two solutions is precise: the derivative of the TAC- 
maker solution has almost linear course while that of the 
|MfcA| solution reveals plantigrade shape. This result allows 
us to assume that our solution of the direct problem for 
the planet transit in the case of linear limb-darkening law is 
more accurate than that of M&A 

For more detailed comparison of the two approaches we 
analyzed the differences (residuals) between the flux values 
of the TAC-maker solution and the IMfcAl solution for the 
same configuration parameters and phases. We established 
that they depend on different parameters (limb-darkening 
coefficients, inclination, phase, planet radius, planet tem- 
perature). 

For linear limb-darkening law the residuals vary with 
the phase by oscillating way around level (Fig. [9| which 
analysis led us to the following conclusions. 

(a) The frequencies and amplitudes of the oscillating resid- 
uals depend on the limb-darkening coefficients (Fig. [9 
The residuals are zero for u = 0.0 that means that the 
solution is precise for uniform star. 



top). 



Mfc~Al 



(b) The frequencies of the residual oscillations decrease to 
the central part of the transit for all orbital inclinations, 
especially for small inclinations (Fig. [9j bottom) . 

(c) The amplitudes of the residuals are bigger for partial 
transits than for total ones (Fig. [9} bottom). 

(d) The residual values increase with the planet radius (as 
it should be expected). 



3.1.2 Comparison for nonlinear limb- darkening laws 

We present comparison with those nonlinear limb-darkening 
laws which are available in the code OCCULTNL. 

(1) The comparison of the synthetic light curves generated 
for the same configuration parameters and quadratic limb- 
darkening law revealed that those produced by the code OC- 
CULTNL exhibited (Fig. 



r.= 1, o = 10r., r pl =0.1r. 
i=90, u=0.2, linear law 
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Figure 8. Comparison of the new solution with that of Mandel 
|fc Agol| p002t for the same parameters: i = 90°, a = 10, R a = 1, 



0.1, T a = 5000 K, T p = 0. Top panel: The coincidence of 



10 1: (i) small-amplitude oscillations 



the two solutions for the linear limb-darkening law and moderate 
precision of the |MfcA| method; Second panel: The difference of 
the two solutions for small part of the transit in big scale; Third 
panel: The first derivatives of the two solutions for small part 
of the transit in big scale; Bottom panel: The residuals for the 
chosen part of the transit in big scale. 



similar to those for the linear limb-darkening law; (ii) the 
fluxes into the transit are systematically smaller than ours 
(the underestimation increases with the coefficient U2 of the 
quadratic term). 

(2) The comparison of the synthetic light curves gener- 
ated for the same configuration parameters and squared root 
limb-darkening law revealed that those produced by the code 
OCCULTNL suffered from the same disadvantages as those for 
quadratic limb-darkening law (Fig. |10[ ). 

The detailed comparison (validation) of our approach 
with the wide-spread M&A method demonstrated that the 
difference between them did not exceed 10 -6 for the con- 
sidered three types limb-darkening laws. This leads to the 
conclusion that for variety of model parameters the new 
approach gives reasonable and expected synthetic transits 
which practically coincide with those of the [Mfc A] solution. 



3.2 Precision vs computational speed 

There is a possibility to increase the precision of the syn- 
thetic transits generated by the code OCCULTNL (E. Agol, 
private communication). Its convergence criterion is the 
maximum change to be less than the product of two mul- 
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Figure 9. Residual oscillations in the case of linear limb- 
darkening law: top for i = 90° and different linear limb-darkening 
coefficients; bottom for different orbital inclinations and u = 1.0 
(all residuals oscillate around level but are shifted vertically for 
a good visibility). 
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Figure 10. Quality of the synthetic transit for different values of 
the precision parameters for the case of quadratic limb darkening 
law and equal coefficients. Top and middle panels: residuals for 
the code OCCULTNL corresponding to different values of PP; bot- 
tom panel: residuals for the code TAC-MAKER corresponding to 
different values of PPtac ■ 



(n is an integer) and the 



tipliers, the factor PP = 10" 
transit depth. 

We established that the decreasing of PP fro m its de- 
fault value 10~ 3 leads to the following effects (Fig. 10 top 
and middle panels): (i) decreasing of the amplitudes of the 
oscillating residuals; (ii) increasing of their frequencies; (iii) 
decreasing of the underestimation of the fluxes for the whole 
transit. 

On the other hand the default value of PPtac is 
1.4 x 10~ 8 (as it is defined for the package SCIPy). Its in- 
creasing leads to some small imperfections of the numerical 
calculations made by our code (Fig. 10 bottom). 



The two considered approaches of the direct problem 
solution lead to different definitions of the precision param- 
eters of the corresponding codes. Therefore, it is not rea- 
sonable to compare the precision of the generated synthetic 
transits for equal values of the parameters PPtac and PP. 
Though, we established that the two methods allow one to 
reach the desired (practically arbitrary) precision of calcu- 
lations by appropriate choice of the corresponding precision 
parameter. Moreover, for the best precisions, the differences 
between the transit fluxes generated with the two codes are 
below lO -10 . 

The computational speed is another important charac- 
teristic of the codes for transit synthesis, especially when 
they are used for the solution of the corresponding inverse 
problem. In principle, this speed depends on two parame- 



ters: time precision (step in phase) and space precision (size 
of the differential element of the integrals). 

The thorough comparison of the computational speeds 
of different codes for the solution of the same problem re- 
quires these codes to be written in the same languages, to 
be run on the same computers and for the equal model pa- 
rameters. That is why it is not reasonable to compare the 
computational speed of our code with the others. 

Nevertheless, we made some tests which revealed the 
following results. 

(a) The reducing of PP (i.e. the increasing of the precision 
of the synthetic transit) by one order leads to 3-5 times 
increase of the computational time while the reducing of 
PP from 10" 3 to 10" 8 (with five orders) requires around 270 
times longer computational time for the code OCCULTNL. 

(b) The decreasing of PPtac from 10" 1 to 10" 8 (with 
seven orders) requires around 15 times longer computational 
time. 

(c) The formal comparison of the absolute computational 
speeds of the Python code TAC-maker and IDL code OC- 
CULTNL revealed that our code is slightly faster for high- 
precision calculations while for the low-precision calculations 
the code occultnl is faster than the code TAC-maker. 
From these results as well as from the consideration that 
each IDL code is faster than its Python version we con- 
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elude that the code TAC-maker possesses good computa- 
tional speed for the high-precision calculations. 



(d) The new subroutine exofast_occultquad (Eastman 
et al. 2013 1 is around 2 orders faster than its progen- 
itor OCCULTQUAD. Moreover, we established that EXO- 
FAST_OCCULTQUAD (with IDL, Fortran and Python ver- 
sions) is considerably improved version of OCCULTQUAD 
not only concerning the computational speed but also con- 
cerning the precision (we had found some bugs of OC- 
CULTQUAD). Although the computational speed of EXO- 
fast_occultquad is higher than those of TAC-maker 
and OCCULTNL, it should be remembered that EXO- 
FAST_OCCULTQUAD can produce transits only for quadratic 
limb-darkening law. 



3.3 A possibility to fit the planet temperature 

The previous transit solutions consider the planet as a black 
lens and the corresponding codes do not fit the planet tem- 
perature T p . As a result only a formal parameter, the equi- 
librium temperature of the planet depending on the stellar 
heating (i.e., on stellar temperature and planet distance), 
can be calculated (out of the procedure of transit modelling). 

The increasing precision of the observations raises the 
problem of the planet temperature, i.e. to study the effect 
of the planet temperature on the transit and to search for 
a possibility to determine this physical parameter from the 
observations. 

The new approach allowed such tests to be directly car- 
ried out because T p is an input parameter of its own (un- 
like the previous methods). As a result, we established that 
increasing of the planet temperature from K to 1000 K 
causes the decreasing of the transit depth up to 10 -6 while 
the increasing of the planet temperature from K to 2000 K 
leads to shallower transit up to 10 -5 . The contribution of the 
planet temperature on the transit depth increases both with 
the ratio R p /R s and limb-darkening coefficients. Moreover, 
it rapidly increases with wavelength. 

The obtained estimations revealed that the effect of the 
planet temperature is lower than the precision of the present 
optical photometric observations, even than that of the Ke- 
pler mission. Hence, the determination of the real planet 
temperature from the observed transit in optics is postponed 
for the near future. 

However, the sensibility of the observations at longer 
wavelengths to the planet emission is higher. Recently, there 
have been reports that the observed IR fluxes of some hot 
Jupiter systems during occultation are higher than those 



corresponding to their equilibrium temperatures ( Gillon et 
aLl[20Tj9l IGibson et al.||2010| ICroll et al. 112010b . These re- 



sults imply that the determination of the planet temperature 
from the observed transits is forthcoming. Until such obser- 
vational precision is reached the code TAC-maker could be 
used for theoretical investigations of the effects of different 
thermal processes (as Ohmic heating, tidal heating, internal 
energy sources, etc.) on the planet transits. 



are quite complex to allow analytical descriptions. Instead 
of analytical solutions we are now able to use fast numerical 
computations. 

This paper presents a new solution of the direct problem 
of the transiting planets. It is based on the transformation of 
the double integrals describing the light decrease during the 
transit to linear ones. We created the code TAC-maker for 
generation of synthetic transits by numerical calculations of 
the linear integrals. 

The validation of our approach was made by compar- 
ison with the results of the wide-spread |M&A| method for 
modelling of planet transits. It was demonstrated that our 
method gave reasonable and expected synthetic transits for 
linear, quadratic and squared-root limb-darkening laws and 
arbitrary combinations of input parameters. 

The main advantages of our approach for the planet 
transits are: 

(1) It gives a possibility, for the first time, to use an arbi- 
trary limb-darkening law /(iij,/i) of the host star; 

(2) It allows acquisition of the stellar limb-darkening coef- 
ficients from the transit solution and comparison with the 
theoretical values of Claret (2 004[ ); 

(3) It gives a possibility, in principle, to determine the 
planet temperature from the observed transits. Our esti- 
mations reveal that the effect of the non-zero planet tem- 
perature to the transit depth is lower than the precision of 
the present optical photometric observations. However, the 
higher sensibility of the observations at longer wavelengths 
(IR) to the planet emission implies that the determination of 
the planet temperature from the observed transits is forth- 
coming. 

These properties of our approach and the practically ar- 
bitrary precision of the calculations of the code TAC-maker 
reveal that our solution of the planet transit problem is able 
to meet the challenges of the continuously increasing photo- 
metric precision of the ground-based and space observations. 

We plan to build an inverse problem solution for the 
planet transit (for determination of the configuration pa- 
rameters) on the basis of our direct problem solution and 
by using the derived simple analytical expressions in this 
paper to obtain initial values of the fitted parameters. 

The code TAC-maker is available for free download 
from the Astrophysics Source Code Librarjj^Jor its own sit^] 
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4 CONCLUSIONS 

Although the analytical solutions of the problems are very 3 http://asterisk.apod.com/wp/ 

useful, the modern astrophysical objects and configurations 4 http://astro.shu-bg.net/software/TAC-maker 
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